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Abstract. In AI2O3 suspensions, depending on the experimental conditions very different mi- 
crostructures can be found, comprising fluid like suspensions, a repulsive structure, and a clustered 
microstructure. For technical processing in ceramics, the knowledge of the microstructure is of im- 
portance, since it essentially determines the stability of a workpiece to be produced. To enlighten 
this topic, we investigate these suspensions under shear by means of simulations. We observe cluster 
formation on two different length scales: the distance of nearest neighbors and on the length scale 
of the system size. We find that the clustering behavior does not depend on the length scale of 
observation. If inter-particle interactions are not attractive the particles form layers in the shear 
flow. The results are summarized in a stability diagram. 

PACS numbers; 82.70.-y, 47.ll.-j, 02.70.Ns, 77.84.Nh 



I. INTRODUCTION 

Colloid science is a very fascinating research 
field, gaining more and more importance in 
the last years. It closely connects physics, 
chemistry, material science, biology, and several 
branches of engineering technology. According 
to its key role in modern science a considerable 
amount of research has been performed to de- 
scribe colloidal suspensions from a theoretical 
point of view and by simulations [l|, 0, H, 0, H, @| 
as well as to understand t he p article-particle 
interactions 0, ol. 11^ . [ill . Il2j . the phase be- 
havior [1^, [ij, flsTllGll. the relevant processes 
on the microscale and their influence on macro- 
scopic parameters [13, [H, El]- Colloidal sus- 
pensions are in fact complicated systems, since 
different time and length scales are involved. 
The particle sizes are on a mesoscopic length 
scale, i.e., in the range of nanometers up to 
micrometers. In systems of particles sized on 
this length scale Brownian motion often cannot 
be neglected. Depending on the particle sizes, 
materials, and concentrations, different interac- 
tions are of relevance and often several of them 
are in a subtle interplay: electrostatic repulsion, 
depletion forces, van der Waals attraction, hy- 
drodynamic interaction, and Brownian motion 
are the most important influences. The prop- 
erties of the suspension are strongly depend- 
ing on the balance of the microscopic forces be- 
tween the particles. Especially for industrial 
processes, where one needs to optimize certain 
material properties a detailed understanding of 
the relevant influences is needed. The stability 
of different microstructures and especially the 
clustering process are key properties which are 



of interest. 

In our work we investigate these properties, 
focusing on AI2O3 particles of diameter 0.37 fim 
suspended in water. This is a widely used ma- 
terial in ceramics [13, HH . We have developed a 
simulation code for a Brownian suspension [2^ 
and have adjusted the simulation parameters 
so that the simulation corresponds quantita- 
tively to a real suspension such that experimen- 
tal data can be compared directly. The diffu- 
sion coefficient, sedimentation velocity 22], and 
the viscosity of the suspension can be repro- 
duced We also have tested the influence 
of polydispersity and found that its influence 
on the results is small. It is much more im- 
portant to choose the correct mean size of the 
particles [H. For AI2O3 suspensions attractive 
van der Waals forces are important for the be- 
havior of this material. Electrostatic repulsion 
of the charged particles counteracts the attrac- 
tion and can prevent clustering depen ding on 
the particle surface charge. In Ref. |23| we 
have presented how one can relate parameters 
of DLVO potentials [3, [1] with experimental 
conditions. In the experiment one can control 
the pH-value and the salt concentration. The 
latter can be expressed by the ionic strength /, 
which is an effective concentration of all ions 
present in the solution. Both, the pH-value 
and the ionic strength, influence the charge of 
the colloidal particles. We have shown that for 
not too strongly attractive forces one can ob- 
tain reasonable quantitative agreement with ex- 
perimental results: by adjusting the lubrication 
force, which represents the short range hydro- 
dynamics, we are able to reproduce rheological 
data of a real suspension [23| . 
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Three regimes can be identified and plotted 
in a stability diagram (23l |. which we want to 
investigate here in more detail: A clustered 
regime, in which particles aggregate to clusters, 
a fluid-like and stable stable suspension and a 
repulsive region, for which the microstructure 
is similar to the ones known from glassy sys- 
tems. From our previous work we know that 
our model works well, even quantitatively, in 
the suspended regime of the stability diagram 
and close to the borders between the different 
microstructures. In this paper we extend our 
investigations to different pH-values, deeper in 
the clustered regime, and to the repulsive struc- 
ture. We expect to gain insight to the micro- 
scopic structure on a qualitative level. 



II. SIMULATION METHOD 

Our simulation method consists of two parts: a 
Molecular Dynamics (MD) code, which treats 
the colloidal particles, and a Stochastic Rota- 
tion Dynamics (SRD) simulation for the fluid 
solvent. 

In the MD part of our simulation the col- 
loidal particles are represented by monodisperse 
spheres. We include effective electrostatic inter- 
actions and van der Waals attraction, known as 
DLVO potentials 0, [1] i ^ lubrication force and 
Hertzian contact forces. DLVO potentials are 
composed of two terms, the first one being an 
exponentially screened Coulomb potential due 
to the surface charge of the suspended particles 



On these grounds we have explored the sta- 
bility diagram of AI2O3 suspensions. The par- 
ticles are uncharged close to the so called "iso- 
electric point" at pH = 8.7. There, for all ionic 
strengths the particles form clusters. For lower 
pH-values particles can be stabilized in solution, 
because they are charged. For low pH-values, 
low salt concentrations, and high volume frac- 
tions a repulsive structure can be found. 



In the following section we shortly describe 
our simulation method. After that we discuss 
the criteria we apply to characterize the mi- 
crostructures. We utilize the pair correlation 
function and the structure factor to character- 
ize the clustering behavior. Both of them in 
principle contain the same information, but we 
concentrate on certain peaks in either of them. 
Each peak in the correlation function and in the 
structure factor corresponds to a certain length 
scale and we chose either the correlation func- 
tion or the structure factor, depending on which 
of the two quantities is more suitable under nu- 
merical criterions to observe on a given length 
scale. In the section thereafter we describe our 
simulation setup. In the results section we start 
with the discussion of the correlation function 
and the structure factor. Additionally, we eval- 
uate the so-called demixing parameter [23 |. 
To characterize the repulsive region we evaluate 
the mean squared displacement (MSD), which 
shows a plateau, if the particle motion consists 
of different processes acting on well separated 
time scales. Finally, the results are summarized 
in a stability diagram for our Al2 03-suspension. 
It shows the behavior of the suspension in an 
intuitive way and helps to design industrial pro- 
cesses using this material. 



Coul 



x^exp( 



2+Kd 
l + Kd 



MbT tanh ^ 



(1) 

r is 
e is 



where d denotes the particle diameter and 
the distance between the particle centers, 
the elementary charge, T the temperature, fcs 
the Boltzmann constant, and z is the valency 
of the ions of added salt. £0 is the permittivity 
of the vacuum, = 81 the relative dielectric 
constant of the solvent, k is the inverse Debye 
length defined by = SttIbI, with the ionic 
strength / and the Bjerrum length £b = 7 A. 
The first fraction in Eq. ([1]) is a correction to the 
DLVO potential (in the form used in Ref. psj). 
which takes the surface curvature into account 
and is valid for spherical particles 26]. 

The effective surface potential ( is the elec- 
trostatic potential at the border between the 
diffuse layer and the compact layer. Smolu- 
chowski related it to the electrokinetic mobil- 
ity of the particle as /i = C^oSr/?? Hzl- The ( 
potential can be related to the pH-value of the 
solvent with a so-called 2pK charge regulation 
model [13. 

The Coulomb term competes with the second 
part of the DLVO potential which is given by 
the attractive van der Waals interaction 



Vvdw= — if 
+2\n 



(2) 



An = 4.76 • 10"^° J is the Hamaker con- 
stant [2^. The DLVO potentials exhibit two 
minima, one primary minimum, where the par- 
ticles touch each other. Eq. ([2]) diverges where 
in reality the primary minimum of the poten- 
tial is located. Therefore we model the primary 
minimum by cutting off the DLVO potentials 
for small separations and substituting them 
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by a parabola which is connected in continu- 
ously difFerentiable manner to the DLVO po- 
tential. This cutoff is made at distances of sev- 
eral nanometers, where the potential has alredy 
reached negative values, i.e. where the poten- 
tial decends to the primary minimum. The pri- 
mary minimum is separated by a potential bar- 
rier from the secondary minimum, which oc- 
curs at larger particle distances and which is 
less deep. For low salt concentrations, i.e., large 
Debye screening lengths, the electrostatic repul- 
sion overcompensates the van der Waals attrac- 
tion and the secondary minimum disappears. 

Long range hydrodynamic interactions are 
taken into account in the simulation for the 
fluid as described below. This can only repro- 
duce interactions correctly down to a certain 
length scale. On shorter distances, a lubrica- 
tion force has to be introduced explicitly in the 
Molecular Dynamics simulation: 
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with the relative velocity Vi-ei,_L projected on the 
connecting line of the particle centers, t] is the 
dynamic viscosity of the fluid. Our implemen- 
tation of the interaction contains cutoff radii 
to connect it to the long-range hydrodynamics 
and to avoid numerical instability for particles 
touching each other 22, 23,]. 
To avoid that the particles penetrate each 
other, we are using a Hertz force described by 
the potential 



Kid-rf^ if r<d, (4) 



where K is an interaction constant and addi- 
tionally a damping term 



Damp 



(5) 



with a damping constant (3. For the integration 
of the translational motion we utilize a velocity 
Verlet algorithm [1^. 

For the simulation of a fluid solvent, many 
different simulation methods have been pro- 
posed: Stokesian Dynamics (SD) d, [2^ f30l . 
Accelerated Stokesian Dynamics (ASD) [3l|, 
[3^ . pair drag simulations Brownian Dy- 
namics (BD) [25l. [33| . Lattice Boltzmann 
(LB) P, [2, 0, l35l. and Stochastic Rotation Dy- 
namics (SRD) [22. [3^ . These mesoscopic 
fluid simulation methods have in common that 
they make certain approximations to reduce the 
computational effort. Some of them include 
thermal noise intrinsically, or it can be included 
consistently. They scale differently with the 
number of embedded particles and the complex- 
ity of the algorithm differs largely. 



We apply the Stochastic Rotation Dynamics 
method (SRD) introduced by Malevanets and 
Kapral [38|, [3^. It intrinsically contains fluc- 
tuations, is easy to implement, and has been 
shown to be well suitable for simulations of col- 
loidal and polymer suspensions [H, [13, [H, [13, 
[40 . 41, 42] and very recently for star-polymers 
in shear flow [i^. The method is also known 
as "Real-coded Lattice Gas" ]36l] or as "multi- 
particle-colhsion dynamics" (MPCD) [3]. It 
is based on so-called fluid particles with con- 
tinuous positions and velocities. A streaming 
step and an interaction step are performed al- 
ternately. In the streaming step, each particle 
i is moved according to 



Y,{t + T)=v,{t)+T^ri{t), 



(6) 



where Yi{t) denotes the position of the parti- 
cle i at time t and t is the time step. In the 
interaction step the fluid particles are sorted 
into cubic cells of a regular lattice and only the 
particles within the same cell interact among 
each other according to an artificial interac- 
tion. The interaction step is designed to ex- 
change momentum among the particles, but at 
the same time to conserve total energy and to- 
tal momentum within each cell, and to be very 
simple, i.e., computationally cheap: each cell 
J is treated independently. First, the mean ve- 
locity vLjit') — jv^jpj '}lf=i ^ is calculated. 
Njit') is the number of fluid particles contained 
in cell j at time t' = t + t. Then, the veloci- 
ties of each fluid particle in cell j are rotated 
according to 

v,(t + T) =u,(t') + 0,(<')-[v.W-Uj(i')]. (7) 

Oj(t') is a rotation matrix, which is indepen- 
dently chosen at random for each time step and 
each cell. We use rotations about one of the co- 
ordinate axes by an angle ±a, with a fixed. The 
coordinate axis as well as the sign of the rota- 
tion are chosen at random, resulting in 6 pos- 
sible rotation matrices. To remove anomalies 
introduced by the regular grid, one can either 
choose a mean free path of the order of the cell 
size or shift the whole grid by a random vec- 
tor once per SRD time step as proposed by Ihle 
and KroU [H,!!^. 

Three different methods to couple the SRD 
and the MD simulation have been introduced 
in the literature. Inoue et al. proposed a 
way to implement no slip boundary condi- 
tions on the particle surface [sl]. Padding 
and Louis very recently came up with full slip 
boundaries, where the fluid particles interact 
via Lennard-Jones potentials with the colloidal 
particles [13] ■ Falck et al. [i^ have developed a 
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"more coarse grained" method which we use for 
the simulations of the present paper and which 
we descibe shortly in the following. 

To couple the colloidal particles to the fluid, 
the colloidal particles are sorted into the SRD 
cells and their velocities are included in the ro- 
tation step. One has to use the mass of each 
particle — colloidal or fluid particle — as a weight 
factor when calculating the mean velocity 

N,{t') 

with Mj{t') = ^ m„ (9) 

i=l 

where we sum over all colloidal and fluid parti- 
cles in the cell, so that Nj{t') is the total num- 
ber of both particles, fluid plus colloidal ones. 
ruk is the mass of the particle with index i and 
Mj{t') gives the total mass contained in cell 
j at time t' — t + t. In summary, the algo- 
rithm for the fluid simulation is described by 
Eqns. (ini)-®- To some of our simulations we 
apply shear. This is realized by explicitly set- 
ting the mean velocity Uj to the shear velocity 
in the cells close to the border of the system. 
Both, colloidal and fluid particles, are involved 
in this additional step. A thermostat is applied 
to remove the energy introduced to the system 
by the shear force. We have described the sim- 
ulation method in more detail in Refs. [2^. [2^. 

III. BACKGROUND 

In this paper we examine the microstructures 
obtained in our simulations for different con- 
ditions. We vary the pH- value and the ionic 
strength /. The shear rate 7 as an external in- 
fluence is varied as well. We classify the mi- 
crostructures in three categories: suspended, 
clustered, and repulsive. In the suspended case, 
the particles can move freely in the fluid and do 
not form stable clusters. In the clustered regime 
the particles form clusters due to attractive van 
der Waals forces. These clusters can be teared 
apart if shear is applied. In some of our simula- 
tions the clusters are very weakly connected and 
at small shear rates they are not only broken up 
into smaller pieces, but they dissolve to freely 
moving individual particles. In this case, we 
assign the microstructure to the suspended re- 
gion, although in complete absence of the shear 
flow clusters are formed. At the borders be- 
tween the different regimes in fact no sharp 
transitions can be observed. The DLVO forces 
rather steadily increase and compete with the 



hydrodynamic interactions. Accordingly, in ex- 
periments one cannot observe a sudden solidifi- 
cation, but a steadily increasing viscosity when 
leaving the suspended regime i23| . 
Similarly as for attractive forces, repulsive in- 
teractions can restrict the mobility of the par- 
ticles. If this happens, the mean squared dis- 
placement of the particles shows a pronounced 
plateau, as it can be found in glassy systems. 
However, we speak of a "repulsive structure" , 
because the change of the viscosity is not as 
strong as in glasses, where it often changes by 
many oders of magnitude, when the glass tran- 
sition is approached. In addition, to claim a 
system shows a glassy behavior would require 
to investigate the temperature dependence of a 
typical time (e.g. particle diffusion time) and to 
show its divergence as the glass temperature is 
approached. This is difficult to do in the frame- 
work of our simulation model [l^l and therefore 
we prefer to speak about a "repulsive structure" 
which might be identified as a colloidal glass in 
future work. 

In this paper we would like to emphasize the 
analysis of the microstructure for different con- 
ditions. Our aim is to reproduce a so-called sta- 
bility diagram by simulations. The stability di- 
agram depicts the respective microstructure de- 
pending on the pH- value and the ionic strength 
/. We apply different numerical tools to ana- 
lyze the microstructure in our simulations and 
finally arrive at a stability diagram shown in 
Fig. 1101 which summarizes the results which we 
present in the following sections. 

The first tool we apply for that is the pair 
correlation function 




(see Ref. [2^ P-55), where V is the volume, N 
the number of particles and rij the distance of 
two particles i and j, can be used for a first 
characterization of the system. It shows max- 
ima at certain typical particle distances, e.g. 
there is a nearest neighbor peak, and more com- 
plicated structures at larger distances, which 
we have assigned to typical particle configura- 
tions for small distances [12] ■ Here, we use the 
next-neighbor-peak to analyze our data, which 
provides information of the short-range struc- 
ture of the suspension. To observe the long- 
range structure we move on to its complemen- 
tary quantity: the structure factor, defined by 

1 ^ 

^^^^^N ^ exp(jk.ri„0, (11) 

/,m— 1 
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where N is the number of particles, and r;™ 
is the vector from particle I to particle m. i 
denotes the imaginary unit here. The structure 
factor is defined in k-space and it is related to 
the pair correlation function in real space by a 
three dimensional Fourier transform: 



S{k) - 1 



J drexp(ik • r)pg(r), 



(12) 



with the density p. Therefore, in principle the 
structure factor contains the same information 
as the pair correlation function. However, due 
to numerical reasons and our implementation of 
shear boundary conditions it is easier to observe 
the long-range structure in the structure factor 
than in the pair correlation function. For our 
evaluation it is important that in a finite sys- 
tem k-vectors are restricted to discrete vectors 
2^ (^jf^ , ^ , jf-^ with Lx.y^z being the system 

size in a;-, y-, and z-direction and Ux^y.z G 
We evaluate ^(k) for these k-vectors and use 
|S'(k)| for further analysis. Since we do not 
evaluate the anisotropy of the structure factor 
in this study, we average over all possible ori- 
entations of the k-vector. We sort the absolute 
values of the k-vectors into intervals and av- 
erage within these intervals over the values of 
the structure factors. Let us now discuss some 
typical features of the structure factor, corre- 
sponding to typical length scales present in the 
system. 

One typical length for dense systems is about 
one particle diameter, the distance particle cen- 
ters have to keep so that they do not overlap. 
The corresponding peak of the structure factor 
is the nearest neighbor peak at fc = In a 
single crystal the peak would be sharp, since 
the particle positions are well defined, whereas 
in our case of a suspension there is a certain 
disorder, which broadens the peak. 

Similar to the nearest neighbor peak another 
peak can be detected at twice the k-vector, 
which corresponds to a distance of one parti- 
cle radius. This does not necessarily mean that 
there are particles whose centers are only one 
particle radius apart. It only means that for a 
certain k-vector the addends in Eq. ([TT|) do not 
cancel out each other. 

For small k-vectors there is another peak, or 
in our case an increase of the structure factor 
towards low k-vectors. The length scale cor- 
responds to the size of the whole system. If 
large clusters are formed, this can be seen in an 
increase of the low-k-peak, since the contribu- 
tions on this length scale do not cancel out each 
other anymore. 

As already mentioned, later in this paper we 
focus on the low-k-peak which corresponds to 



a length scale of the size of the system. 



IV. SIMULATION SETUP 

In this study the colloidal particles are rep- 
resented by three dimensional spheres oi d = 
0.37 /zm in diameter. This is the mean diam- 
eter of the particles used in the experiments 
to which we refer in Ref. [1^. We have sim- 
ulated a small volume, 24 d = 8.88 /j,m long 
in x-direction, which is the shear direction, 
and 12 d = 4.44 ^m long in y- and z-direction. 
We have varied the volume fraction between 
$ = 10 % (660 particles) and $ = 40 % (2640 
particles). Most of the simulations were per- 
formed at $ = 35% (2310 particles). We use 
periodic boundaries in x- and y-direction and 
closed boundaries in z-direction Shear is 

applied in x-direction by moving small zones of 
particles and fluid close to the wall with a given 
shear velocity. The xy-plane is our shear plane. 
For simulations without shear, to achieve the 
best comparability, we use the same boundary 
conditions and just set the shear rate to 7 = 0. 
In addition we have performed simulations with 
two different shear rates: with 7 = 100/s and 
with 7 — 500/s. For better comparability to 
other works the Peclet number 



Pe = 67n]R^/kBT, 



(13) 



which expresses the importance of Brownian 
motion with respect to the shear flow, is useful. 
The Peclet numbers in our shear simulations 
are Pe ~ 3 and Pe = 15, respectively. 



V. RESULTS AND DISCUSSION 

First, we focus on simulations without shear, 
where one can predict intuitively, what should 
happen. Qualitatively the results are similar to 
our earlier work |22l | , but the quantitative rela- 
tion between the pH-value and the potentials is 
new. The relation was presented in Ref. psj . 
but here we apply it to different cases and we 
focus more on the characterization of the mi- 
crostructure. However, given the particle par- 
ticle interaction potentials, the microstructure 
in equilibrium can be predicted easily, at least 
on a qualitative level. But, the matter changes 
and gets more sophisticated, when shear is ap- 
plied and an interplay between shear flow and 
particle particle interactions becomes responsi- 
ble for the resulting microstructure. 
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A. Correlation function 

For constant ionic strength / = Smmol/l the 
local microstructure can be examined using 
the correlation function. Depending on the 
pH-value the behavior of the system changes 
from a repulsive structure around pK = 4 to 
a stable suspension around pH = 6 towards a 
clustered region if the pH-value is further in- 
creased, until the isoelectric point is reached at 
pH — 8.7. There clustering occurs in any case, 
independent on the ionic strength. This can be 
seen in the structure of the correlation function: 
electrostatic repulsion prevents clustering (at 
pH = 4). Particles are suspended, and there is 
no fixed long range ordering in the system. The 
correlation function (Fig.[T]) shows a maximum 
at a typical nearest neighbor distance slightly 
above ^ = 1 with d denoting the particle diam- 
eter, then in the layer of next neighbors small 
correlations can be found (at ^ = 2). For larger 
distances the correlation function is rather con- 
stant. 

When the pH-value is increased, the surface 
charge is lower, which at first causes the parti- 
cles to approach each other more closely. The 
maximum of the correlation function is shifted 
to smaller distances (see Fig.[Tl note that the 
curves are shifted vertically in the plot by a 
factor of 3 for better visibility). Then, van 
der Waals attraction becomes more important 
and clustering begins. One can see this in the 
correlation function where a sharp structure at 
particle distances between 1.5 and 2 particle di- 
ameters occurs. In a solid like cluster the posi- 
tion of the next neighbor is fixed more sharply 
than in the suspension, consequently the near- 
est neighbor peak becomes sharper, too, and 
its height is increased. Close to the isoelec- 
tric point (pH = 8.7) the barrier between pri- 
mary and secondary minimum disappears. The 
particles, once clustered, cannot rearrange any- 
more, and therefore the correlations to the next 
neighbors become less sharp again (compare the 
cases of pH = 8.7 and pH = 7.7 in Fig.^at the 
positions denoted by the arrows). 

Instead of varying the pH- value, one can also 
vary the ionic strength to achieve similar ef- 
fects. Increasing the ionic strength, experi- 
mentally speaking "adding salt" decreases the 
screening length 1/k and therefore the attrac- 
tive forces become more important: the parti- 
cles start to form clusters. On the contrary, if 
the ionic strength is decreased-experimentally 
speaking a dialysis step is performed-the elec- 
trostatic repulsion prevents cluster formation 
and, for sufficiently strong repulsion, long range 
correlations occur as soon as the range of the 
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Figure 1: Dependence of the particle correlation 
function on the pH value, / = 3mmol, 7 = 0/s 
$ — 35%. The plots for four different pH- values 
are shifted against each other for better visibility 
by a factor of 3. For pH — 4 the particles are not 
clustered. Hence the structure at ^ = 2 is less 
sharp than in the other three curves of the plot 
and the nearest neighbor peak (at ^ = 1) is broad. 
For pH — 6.5 slight clustering starts, the structures 
become sharper. For pH = 7.7 strong cluster for- 
mation is reflected in very sharp structures. For 
pH = 8.5 electrostatic repulsion nearly disappears 
so that no barrier between primary and secondary 
minimum exists anymore. The particles cannot re- 
arrange anymore, and therefore the structures la- 
beled by the arrows become smoothened compared 
to the case of pH = 7.7. 



repulsion reaches the mean particle distance. 

In Fig.[2]one can see clustering in the primary 
minimum of the potential as well as in the very 
shallow secondary minimum. The correlation 
function plotted there is evaluated after 1000 
SRD time steps in a simulation for pH — 7 and 
/ = 3 mmol/1 at a volume fraction of $ = 35% 
sheared with 7 = 100/s. Note that the shear 
flow is not of essential importance here, but it 
supports the particles to overcome the potential 
barrier between the primary and the secondary 
minimum. The dotted lines in Fig.[2]denote the 
zero line. 

A secondary minimum only exists if the 
screening length of the electrostatic repulsion is 
short enough, i.e., if the ionic strength is large 
enough. According to the depth of the poten- 
tial, clustering in the primary minimum is as- 
sociated with much stronger forces than in the 
secondary minimum. 

The effects described up to here can be ob- 
served with or without shear qualitatively in 
an analogous manner. If the suspension is 
sheared clustering occurs at higher pH-values 
and the peaks found in the correlation func- 
tion are slightly broadened, because the relative 
particle positions are less fixed. But a new fea- 
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Figure 2: Correlation function, 7 = 100/s, pH = 7, 
7 = 3mmol/l, $ = 35%: one finds clustering in the 
primary and the secondary minimum. Note: The 
vertical axis is logarithmic, and a constant offset 
has been added before plotting. Thereby, the co- 
incidence of the minima in the potential with the 
maxima in the correlation function becomes more 
apparent. The dotted lines denote the position of 
the base line before shifting. 



ture appears, if a stable suspension of not too 
high volume fraction is sheared. Induced by 
the shear particles arrange themselves in lay- 
ers. Regular nearest neighbor distances in the 
shear plane cause the correlation function to be- 
come more structured even for large distances 
(see Fig. 131). The long range structure of the 
pair correlation function appears after a tran- 
sient time the particles need to arrange them- 
selves in the layered structure. The curves in 
Fig. [3] correspond to the same simulation, but 
for increasing time from bottom to top. The 
arrows from the bottom indicate the position 
of the next-neighbor peak and the next but 
one neighbors (Since the particles do not touch 
each other, the peaks are located not exactly 
at 2 — 1:2,3,...). Within one layer which 
moves as a whole, a hexagonal particle order 
appears, which can be seen in the occurrence of 
a peak at ^ = VS times the position of the near- 
est neighbor peak, i.e., the next neighbor peak 
splits, as indicated by the arrows from above. 
The same applies to the next but one neigh- 
bor peak marked by the second set of arrows 
at 3 < r/d < 4. Shear induced layer formation 
has been found in both, experiments [49l Isoj 
and simulations [5l|, [H, 

We have integrated over the nearest neigh- 
bor peaks, both, the peaks of the primary and 
the secondary minimum, and plotted the inte- 
gral versus pH-value in Fig.m We have cho- 
sen / = 3mmol/l and <i> = 35% and three 
different shear rates: 7 = 0,100 and 500/s. 
We have integrated the correlation function for 



Figure 3: Correlation function for 7 = 500/s, 
pH = 6, / = 0.3mmol/l, $ = 32%: depending 
on the simulation time the peaks indicated by the 
arrows from below split into two (indicated by the 
arrows from above), and long range correlations oc- 
cur (for ^ > 4). This reflects the process of layer 
formation and appearing of a regular (hexagonal) 
order in the layers. From bottom to top several 
states for increasing time are shown. The plots are 
shifted vertically by a factor of 2 for better visibil- 
ity. 



r < 1.215 d, where for all pH- values the po- 
tential in the secondary minimum has a value 
of — i/cBT- In other words, we have captured 
the primary and the secondary minimum of the 
potential for this plot. For low pH- values clus- 
tering (in the secondary minimum) is only pos- 
sible for low shear rates. For high shear rates, 
the hydrodynamic forces do not allow the for- 
mation of stable clusters. For rising pH-values 
the clustering increases, first for the un-sheared 
suspension, at higher pH-values for low shear 
rates (7 = 100/s) and finally for high shear 
rates (7 = 500/s). Remarkably, for pH > 7.5 
the curve for 7 = 100/s shows stronger clus- 
ter formation than the other ones. Particles 
are brought together by the shear flow, so that 
compared to the case of no shear, the clustering 
process is supported here. On the other hand, 
the shear stress may not be too strong, because 
otherwise the clustering process is limited by 
the shear flow again (for 7 ~ 500/s the cluster- 
ing is less pronounced than for 7 = 100/s). 



B. Structure Factor 

The pair correlation function can be used to 
characterize the local order of the microstruc- 
ture on the length scale of the particle size. 
However, to do the characterization on the 
length scale of the system size, we use the struc- 
ture factor S'(k). In principle we could integrate 
the correlation function g{r) over an interval for 
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Figure 4: Nearest neighbor peak (primary and sec- 
ondary minimum of the potential) of the correla- 
tion function / = 3mmol/l, "1> = 35%: For low 
pH-values clustering is prevented by the electro- 
static repulsion. For high pH-values the particles 
form clusters, which is reflected by an increased 
nearest neighbor peak. First, shear prevents clus- 
tering, then depending on the shear rate, cluster 
formation takes place. Low shear rates even sup- 
port cluster formation at high pif-values. 



large r. However, due to our implementation 
of shear (see Ref. [23j), we have to close the 
boundary in z-direction. Therefore we already 
find restrictions for r around half the extension 
of the system in this direction. We find that it 
is easier to handle these finite size restrictions 
by moving on to the structure factor. There, 
we always find a low-k-peak, even if no cluster 
formation takes place in the simulation, fn con- 
trast to experiments, where this peak only ap- 
pears for clustered samples, it reflects the pres- 
ence of a typical length of the system size in the 
simulation. Namely, the finite size of the sim- 
ulation volume is the typical length, which ap- 
pears in the structure factor by the low-k-peak. 
This only produces a constant offset when in- 
tegrating over the low-k-peak as we are going 
to do below in this section. Thus, it is easy to 
handle the influence of the finite system size in 
k-space. 

In Fig. [5] we have plotted several typical 
structure factors of our simulations. For these 
plots the pH- value is fixed to pH = 6. The cases 
a) and b) are sheared with 7 — 500/s at an ionic 
strength of / = 0.3mmol/l. In case a) the vol- 
ume fraction $ = 20% is relatively low. There- 
fore the particles can arrange themselves in lay- 
ers parallel to the shear plane, which move rel- 
atively independently in the shear flow. They 
have a certain distance flxed in space and time. 
This can be seen in a sharp peak at a dimension- 
less k- vector of k = 5.2, which corresponds to a 
distance of 1.2 particle diameters. In fact, this 
is exactly the distance between two neighboring 
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Figure 5: Structure factor for some selected exam- 
ples, with pH = 6 fixed for all plots: 7 = 500/s, 
/ = 0.3mmol/l: a) <1> = 20% and b) $ = 35% , 
7 = 0, / = 25 mmol/1: c) $ = 40% and d) $ = 10% 
. The curves are shifted vertically for better visi- 
bility. In case a) ten layers can be identified in 
the system, resulting in the strong peak close to 
5. But, since the particles in the layers can still 
move freely, there is no 2"'''-order-peak. In case b) 
layers are formed, but particles are moving from 
one layer to the other, disturbing the flow. As a 
result the nearest neighbor peak is much broader. 
Due to the structure in the layers, a 2"''-order-peak 
appears. In case c) the interaction is strongly at- 
tractive, hence the particles approach each other 
and the nearest neighbor peak is shifted to higher 
k-vectors. In case d) the volume fraction is much 
less. The slope of the low-k-peak is much flatter, 
which depicts that the cluster is fractal. 





a) 



b) 





d) 



Figure 6: (Color online) Snapshots of the systems 
analyzed in Fig.[Sl In case a) one can see the layers 
resulting in the sharp peak in the structure factor. 
In case b) the layers are packed closer due to the 
higher volume fraction. Collisions between parti- 
cles of neighboring layers happen more frequently. 
In case c) one big cluster is formed. The particles 
are packed densely. In case d) the fractal nature of 
the system can be seen directly. 
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layers, as one can easily verify by counting the 
layers in a snapshot of the system (Fig.[6K))- 
The particles in the layers do not have a fixed 
distance and therefore no 2'"'-order-peak can 
be observed. 

For case b) the volume fraction is increased to 
$ = 35%. The particle layers are packed more 
densely and therefore the interactions between 
one layer and the neighboring one become rel- 
evant. Particles jump from one layer to the 
other, which disturbs the flow and therefore the 
distance between the layers is not fixed any- 
more. The sharp peak on top of the nearest 
neighbor peak disappears. Instead of that, in 
each layer a regular hexagonal order appears 
and therefore the 2"''-order-peak is much more 
pronounced. 

In case c) the ionic strength is increased to 
/ — 25mmol/l. The inter-particle potentials 
are attractive enough that aggregation takes 
place. In this simulation we did not apply 
shear, therefore one finds only one big cluster 
in the system (compare Fig.|BJ:)). In the clus- 
ter the particles are packed more densely and 
consistently the nearest neighbor peak in the 
structure factor is shifted to larger k-vectors. 
The volume fraction is = 40% in this case. 
In case d) the volume fraction is decreased to 
$ = 10%. The particles still form clusters, 
but their mobility is not high enough to create 
one compact cluster. The system has a frac- 
tal structure (see Fig.[6ji)). This can be seen in 
the structure factor as well: The slope of the 
low-k-peak is flatter in this case compared to 
cases a)-c) . A flatter slope of the low-k-peak is 
typical for structure factors of fractal objects. 
The fractal dimension of the cluster extracted 
from the slope of the low-k-peak is 2.5. In ex- 
periments this relation is often used to deter- 
mine the fractal dimension of a sample: Lat- 
tuada et al. [131 have evaluated the fractal di- 
mension of agglomerates of latex particles from 
the slope of the structure factor. McCarthy 
et al. [55| give an introduction to scattering in- 
tensities at fractal objects, without mentioning 
the structure factor, but their arguments refer 
to the contribution of the structure factor on 
the scattering intensity. The underlying mech- 
anism which is responsible for these structures 
is cluster cluster aggregation [s^. 
In Fig. [7] we show the dependence of the 
low-k-peak of the structure factor on the 
pH-value. Here we have integrated over dimen- 
sionless k-vectors smaller than 3 which means, 
we have captured structures larger than twice 
a particle diameter. A large integral over the 
low-k-peak is due to a large inhomogeneity in 
the system. In one part of the system parti- 
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Figure 7: low-k-peak for different pH-values and 
different shear rates. The ionic strength / is kept 
constant at / = 3 mmol/1 and the volume fraction is 
always $ — 35%. For 7 = 0/s the particles tend to 
cluster in the secondary minimum of the potential. 
This clustering can easily be broken up, if shear is 
applied. If the pH- value is increased, shear cannot 
prevent cluster formation anymore. At low shear 
rates (7 — 100/s) clustering is even enhanced, since 
the particles are brought closer to each other by 
the shear flow. Note that the offset of the plots 
reflects the finite size of our system combined with 
the closed boundary conditions (see the text). 



cles are present and in the other part not. In 
other words, we observe the process of clus- 
ter formation on a length scale of the system 
size. Without shear, particles cluster in the 
secondary minimum for all pH- values. If the 
system is slightly sheared (7 — 100/s) cluster- 
ing is suppressed for low pH-values. Starting 
at pK = 6 cluster formation starts and is even 
supported by the shear flow for pH- values larger 
than 7.5. For large shear rates (7 = 500/s) clus- 
ter formation is suppressed by the shear flow. 
By analyzing the low-k-peak of the structure 
factor one observes on the length scale of the 
system size. The same behavior of the system 
can be seen by analyzing the pair correlation 
function, as we have already shown in Fig.[31 
In that case one analyzes the number of nearest 
neighbors, that means, one observes the length 
scale of a particle diameter. Nevertheless, both 
graphs show the same behavior of the system, 
i.e., we have a consistent picture of the clus- 
ter formation process on the length scale of the 
nearest neighbors and on the length scale of the 
system size. 

Thus we have confirmed that the cluster for- 
mation process is not limited to length scales 
smaller than our system size. This is reflected 
especially by the transition between pK = 7—8 
and its shear rate dependence in the plots in 
Fig. [7] and Fig. 21 There is a strong similarity of 
the two plots, which are obtained by two evalu- 
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ation methods referring to two different length 
scales. This confirms that the plots do not only 
reflect how clusters are formed on the respective 
length scale, but that the clustering process is 
a phenomenon which can be observed on any 
length scale by applying a suitable method to 
characterize it. 



C. Density Inhomogeneity 

Another way to observe the cluster formation 
process is provided by the "demixing parame- 
ter" defined in Ref. [24| as follows: the system 
is divided into n!^ cubes and the particle density 
Pk is evaluated in each cube. Then the demix- 
ing parameter is the mean squared deviation of 
the density 

^r.^Y.^pU--p)\ (14) 
fc=0 

where p is the mean density. For our elongated 
system we modify this definition and use 2?! 
cubes in x-direction, resulting in 271^ cubes in 
total. The demixing parameter implicitly con- 
tains information about the density of the clus- 
ters. It is the density fluctuation of the whole 
system. If the clusters are more compact, voids 
have to appear and the demixing parameter in- 
creases, since the distribution of the local den- 
sity is broadened thereby. 

In Fig. [5] we have plotted the demixing pa- 
rameter \['4 versus time for three different shear 
rates and four different ionic strengths in each 
plot. The pH- value is kept constant at pH = 6 
and the volume fraction at $ = 35%. Without 
shear ^'4 is nearly constant. Even if a cluster 
is formed, it does not move, so that the aver- 
age density in the boxes does not change much. 
For lower volume fractions the time dependence 
without shear is stronger. If shear is applied 
the clusters are deformed so that stronger in- 
homogeneities appear. Depending on the ionic 
strength the demixing parameter increases with 
cluster formation (high ionic strength) or still 
stays constant when the interactions are repul- 
sive (low ionic strength). For 7 = 100/s the 
strongest cluster formation is achieved for the 
highest ionic strength / = 20mmol/l. If the 
ionic strength is decreased, the clustering effect 
decreases as well, until it disappears completely 
at / = 1 mmol/1, where the repulsive regime 
is reached. To see the shear rate dependence, 
compare the plots for 7 = 100/s and 500/s. If 
the shear rate is increased, the cluster forma- 
tion starts faster because the shear flow brings 



the particles faster in contact, but the result- 
ing inhomogeneity of the final state is less than 
for lower shear rates. This effect can be seen 
in Fig.U] and Fig. [71 too. There, also the in- 
tensity of the respective peaks for 7 = 500/s 
is less than for 7 = 100/s. This shows again 
that large shear rates can inhibit cluster for- 
mation, whereas moderate shear rates can sup- 
port the clustering process (the plots in Fig.[8]b) 
for 7 = 100/s are steeper than the ones in 
Fig.[5]a) for 7 = 0/s.) In addition to the infor- 
mation already contained in other quantities we 
have presented, Fig.[8]shows the time evolution. 
One can see three regimes of time evolution: 
for very short times all plots arc nearly con- 
stant. This means that the cluster formation 
has not yet started. Then, the plots for attrac- 
tive potentials raise, which shows the onset of 
cluster formation. For large times the slope of 
the plots decreases, which reflects that the clus- 
ters have been formed and no more "demixing" 
takes place. Without shear (see Fig.[H]a) ), it is 
remarkable that for the largest ionic strength 
the growth of the demixing parameter first be- 
comes constant after 0.3 s, whereas for the other 
cases where clustering occurs, it grows further. 
This reflects that in the case of / = 20 mmol/1 
the attraction is that strong that no reordering 
of the particles is possible anymore. 



D. Repulsive Regime 

To characterize the repulsive regime, we eval- 
uate the mean squared displacement for the 
particles. In Fig. [3] we plot the mean squared 
displacement for different ionic strengths. The 
pH- value is kept constant at pH — 6 and the vol- 
ume fraction is $ = 35% for this plot. Three 
different regimes can be identified. For very 
short times, the ballistic regime: particles move 
on short distances without a notable influence 
by their neighbors. The distances are in the 
order of some percent of the particle diame- 
ter and the times are a few SRD time steps. 
For larger times the particles interact with their 
neighbors and therefore their mobility is limited 
due to collisions with the neighbors. This is 
reflected in the mean squared displacement by 
a plateau of reduced slope, which is the more 
pronounced the more the mobility of the par- 
ticles is restricted. For even larger time scales 
collective motion starts, i.e., clusters or groups 
of particles move, or single particles can escape 
from a cage formed by its neighbors. Depend- 
ing on the ionic strength different effects are 
important and thus the shape of the curve is 
different. For large ionic strengths the particles 
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Figure 8: Demixing parameter for constant pH = 6 and volume fraction <1> = 35 % for three different 
shear rates: -y — 0/s (left), 100/s (center) and 500/s (right). Without shear only very weak demixing 
takes place. For 7 = 100/s the strongest demixing can be observed for the system with the highest ionic 
strength I = 20 mmol/l. For lower ionic strengths the effect decreases until it disappears completely at 
I = 1 mmol/l, where the repulsive regime is reached. If the shear rate is further increased (7 — 50Q/s), the 
clusters are less stable so that the system becomes less inhomogeneous. For / — 5 mmol/l the demixing is 
even suppressed completely. 
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Figure 9: Mean squared displacement at pH = 6 
for different ionic strengths, without shear. One 
can see a ballistic regime for short times, a central 
plateau and a collective long time movement which 
can be a movement of a whole cluster or cage es- 
cape events of single particles. Depending on the 
ionic strength, the central plateau is more or less 
pronounced. A comparison of the plateau for dif- 
ferent simulations can be used to decide, if a certain 
state belongs to the repulsive region of the stability 
diagram. A state well in the liquid microstructure 
should be used as a reference for the comparison. 



form clusters and these clusters may drift or ro- 
tate in the system. Then the collective motion 
is more dominant and the mean squared dis- 
placement grows faster than in single particle 
diffusion. The mean squared displacement does 
not show a plateau, then. But in the repulsive 
regime, the neighbors limit the motion of the 
particles, and the slope of the plateau is flatter, 
i.e., the plateau is even more pronounced, com- 
pared to the suspended case. In the repulsive 
regime the particles tend to arrange themselves 
in layers when shear is applied |23j and long 
range correlations can be found in un-sheared 
systems [22| . 



E. Stability Diagram 

The results of the investigations presented in 
this paper can be summarized in a stability dia- 
gram for our Al203-suspension (Fig.[Tni). Three 
different microstructures can be identified: a 
repulsive structure, a suspended region and a 
clustered region. In contrast to our previous 
work p^ . [23j . we have explored the param- 
eter space more in the repulsive regime and 
deeper in the clustered region. We use the mean 
squared displacement, the demixing parameter 
^ 2JI, the correlation function, and the struc- 
ture factor, to decide to which of the three mi- 
crostructures a certain point in the stability di- 
agram belongs. However, the borders between 
the regions are not sharp and they depend on 
the shear rate. We have indicated the crossover 
regions by the shaded patterns in the stability 
diagram. If the volume fraction is decreased, 
the region of the repulsive structure becomes 
smaller. 

To decide if a state is in the suspended re- 
gion or in the repulsive one of the stability dia- 
gram, we have compared the plots of the mean 
squared displacement for the simulations with- 
out shear. If the plateau was pronounced there, 
we have counted the state among the repulsive 
regime. As a second criterion one can compare 
the pair correlation function. If there are long 
range correlations even though the system is 
not sheared, then the microstructure is the re- 
pulsive one. Finally, the shear force can be used 
to localize the border to the repulsive regime. 
For a given shear rate and a fixed volume frac- 
tion, the shear force depends on the particle 
interactions. If the shear force increases com- 
pared to a state well in the suspended regime, 
the motion of the particles is blocked by the 
electrostatic interaction in the repulsive regime. 

As we have mentioned in Sec. lIIIi without 
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shear weak clustering can be seen in the sus- 
pended case as weh, since there is no barrier 
for the particles to enter the secondary mini- 
mum of the DLVO potential, but the clusters 
can be broken up again very easily. Thus, to 
decide, if a state belongs to the clustered or to 
the suspended regime, we first study the snap- 
shots of the system. If we see no clusters there, 
the clustered regime can be excluded. But, if 
we see clusters, we next check the density of 
the clusters and the onset time of the increase 
of the demixing parameter which is a mea- 
sure for the time it takes to form clusters. Both, 
the density and the time are indications for the 
stability of the clusters. If they grow slowly 
and their density is low, we count the state 
to the suspended regime. The demixing pa- 
rameter shown in Fig. [8] implicitly contains 
information about the density of the clusters. 
It is the density fluctuation of the whole sys- 
tem. If the clusters are more compact, voids 
have to appear and the demixing parameter in- 
creases, since the distribution of the local den- 
sity is broadened thereby. The stability dia- 
gram obtained by these criteria is consistent 
with the results of the simulations with shear 
flow, shown in Fig. 2] and Fig. [71 Especially, the 
increased cluster formation for / = 3mmol/l 
starting between pH = 7 — 8 is reflected in an 
increased nearest neighbor peak in Fig.[4l and 
low-k-peak in Fig.[7| respectively, and in a bor- 
der between suspended and clustered regime in 
Fig. 1101 The repulsive structure for pH = 4 and 
/ — 3mmol/l can not be recognized in Fig. [3] 
and Fig. [71 but in a pronounced layer forma- 
tion. 



VI. SUMMARY AND OUTLOOK 

We have simulated colloidal particles in shear 
flow and investigated how the clustering pro- 
cess due to attractive DLVO potentials is af- 
fected by the hydrodynamic forces. We find a 
consistent behavior on different length scales. 
The nearest neighbor peak of the pair corre- 
lation function has been used to observe the 
direct neighborhood of the particles and the 
low-k-peak of the structure factor to keep track 
of the length scales up to the system size. In 
both cases a suppression of the cluster for- 
mation by the shear flow can be seen at low 
pH-values. For large pH-values low shear rates 
even support the clustering process. In con- 
trast, for high shear rates it suppresses the clus- 
ter formation. We have evaluated the mean 
squared displacement and the demixing param- 
eter [131 in order to draw the stability dia- 
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Figure 10: Stability diagram (plotted for <1> = 35% 
and without shear): depicting three regions: a 
clustered region (filled circles) , a suspended regime 
(open squares), and a repulsive structure (filled 
squares). In the clustered region particles aggre- 
gate which leads to inhomogeneity in the system. In 
the suspended regime, the particles are distributed 
homogeneously in the system and they can move 
freely. In the repulsive regime the mobility of the 
particles is restricted by electrostatic repulsion ex- 
erted by their neighbors. As a result they arrange 
in a local order which maximizes nearest neighbor 
distances. The borders between the regimes are not 
sharp. They depend on the shear rate and on the 
volume fraction. Therefore we have indicated the 
crossover regions by the shaded patterns. The lines 
are guides to the eye. 



gram as given in Fig.lTUl To our knowledge this 
stability diagram for AI2O3 suspensions is re- 
produced quantitatively for the first time from 
simulations. It helps to predict the behavior of 
a real suspension. Our findings on the cluster 
formation process suggest that soft stirring can 
enhance the cluster formation in industrial pro- 
cessing of this material. Further investigations 
can be carried out on the fractal dimension and 
its dependence on the experimental conditions. 
The low-k-peak of the structure factor can be 
used for that. The cluster size distribution 
could as well deliver interesting insights use- 
ful to design industrial processes. Apart from 
that one can apply our algorithm to different 
materials. To do so, one has to change the in- 
teraction constants and especially for the elec- 
trostatic repulsion the calibration we have pre- 
sented in Ref. [1^ are necessary. We expect 
that the stability diagram does not change qual- 
itatively, but the position of the borders will be 
different. 
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